
% Surendra Adhikari
% Jet Propulsion Lab, Caltech 
% adhikari@jpl.nasa.gov

% June 28, 2020. 

% This Matlab code is a companion to our The Cryosphere submission. 
% https://doi.org/10.5194/tc-2020-23 

% Read 'Adhikari2020TCFigure3_README.txt' before you read this file. 

% Here we guide the user to use the BASIC data provided in the following file and compute DERIVED fields, which are supplied in 'Adhikari2020TCFigure3_DERIVED.mat' file. 

% load the BASIC data. 
load('Adhikari2020TCFigure3_BASIC.mat'); 
% A quick reminder of the key variables: 
% BASE => Base of the ice.
% BED => Bedrock elevation. 
% SURFACE => Surface of the ice. 
% TIMES => Simulation times: From AD 2000 to AD 2490. 
% X => X-coordinate 
% Y => Y-coordinate 

% define some constants.
rho_ice = 917;		% kg/m^3 
rho_ocean = 1028; % kg/m^3 
rho_water = 1000; % kg/m^3 

% compute ice thickness 
THICK = SURFACE - BASE; 

% Change in ice thickness relative to AD 2000. 
DEL_H = THICK - THICK(:,1); 

% Level set function. In the present context, it is sufficent to consider equations 1 and 2...
% ...(rather than equations 3 and 4) to define the ocean and land domains, and their interfaces. 
F = THICK + rho_ocean/rho_ice*BED; % Eq. 1 of the paper. 

% Grounded ice masks. Equivalent to Equation 5.  
GROUND_MASK = F; 
GROUND_MASK(GROUND_MASK<0) = 0; 
GROUND_MASK(GROUND_MASK>0) = 1;	

% Floatation height and the height above flotation (HAF) 
H0 = -rho_ocean/rho_ice * min(BED,0);	% Equation 7. 
HF = GROUND_MASK.*(THICK - H0);			% Equation 8.  

% The HAF method. 
% Change in HAF relative to AD 2000 
DEL_HF = HF - HF(:,1);	% Equation 9. 

% New method of computing sea-level contribution from ice sheets. 
% t = AD 2000
% t + dt = future times 
GG = GROUND_MASK(:,1).*GROUND_MASK; 
DEL_HM = DEL_H.*GG + DEL_HF.*(1-GG);	% Equation 11. 
DEL_HV = (1-rho_water/rho_ocean)*(DEL_H - DEL_HF).*(1-GG);	% Equation 12. 
DEL_HS = DEL_HM + DEL_HV;	% Equation 10. 

